3.8.2.1 F_COL

The maximum number of segregating sites per locus in the data set is 7, there are 30) F_COL run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 30 individuals (60 observations).

# 2*number of individuals in data set
N=60

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind60_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind60")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_COL"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all <- list()

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.15772, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.2 F_MIL

The maximum number of segregating sites per locus in the data set is 7, there are 20) F_MIL run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 20 individuals (40 observations).

# 2*number of individuals in data set
N=40

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind40_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind40")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_MIL"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.21342, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.3 F_DER

The maximum number of segregating sites per locus in the data set is 7, there are 15) F_DER run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 15 individuals (30 observations).

# 2*number of individuals in data set
N=30

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind30_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind30")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_DER"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.12052, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.4 F_BUT

The maximum number of segregating sites per locus in the data set is 7, there are 21) F_BUT run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 21 individuals (42 observations).

# 2*number of individuals in data set
N=42

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind42_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind60")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_BUT"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.20954, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.5 F_FRH

The maximum number of segregating sites per locus in the data set is 7, there are 27) F_FRH run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 27 individuals (54 observations).

# 2*number of individuals in data set
N=54

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind54_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind54")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_FRH"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.10915, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.6 F_NIM

The maximum number of segregating sites per locus in the data set is 7, there are 30) F_NIM run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 30 individuals (60 observations).

# 2*number of individuals in data set
N=60

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind60_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind60")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_NIM"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.16563, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.7 F_MKH

The maximum number of segregating sites per locus in the data set is 7, there are 28) F_MKH run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 28 individuals (56 observations).

# 2*number of individuals in data set
N=56

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind56_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind56")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_MKH"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.16758, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.8 F_STN

The maximum number of segregating sites per locus in the data set is 7, there are 23) F_STN run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 23 individuals (46 observations).

# 2*number of individuals in data set
N=46

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind46_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind46")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_STN"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.19666, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.9 F_TOU

The maximum number of segregating sites per locus in the data set is 7, there are 30) F_TOU run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 30 individuals (60 observations).

# 2*number of individuals in data set
N=60

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind60_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind60")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_TOU"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.1764, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.10 F_MRH

The maximum number of segregating sites per locus in the data set is 7, there are 15) F_MRH run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 15 individuals (30 observations).

# 2*number of individuals in data set
N=30

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind30_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind30")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_MRH"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.19097, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.11 F_MER

The maximum number of segregating sites per locus in the data set is 7, there are 31) F_MER run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 31 individuals (62 observations).

# 2*number of individuals in data set
N=62

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind62_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind62")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "F_MER"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.16428, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.12 L_USR

The maximum number of segregating sites per locus in the data set is 7, there are 21) L_USR run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 21 individuals (42 observations).

# 2*number of individuals in data set
N=42

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind42_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind42")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "L_USR"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.12979, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.13 W_USR

The maximum number of segregating sites per locus in the data set is 7, there are 26) W_USR run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 26 individuals (52 observations).

# 2*number of individuals in data set
N=52

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind52_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind42")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "W_USR"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.21578, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.14 S_MIL

The maximum number of segregating sites per locus in the data set is 7, there are 16) S_MIL run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 16 individuals (32 observations).

# 2*number of individuals in data set
N=32

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind32_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind32")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "S_MIL"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.16933, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.15 S_DER

The maximum number of segregating sites per locus in the data set is 7, there are 27) S_DER run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 27 individuals (54 observations).

# 2*number of individuals in data set
N=54

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind54_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind54")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "S_DER"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.14984, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.16 S_BUT

The maximum number of segregating sites per locus in the data set is 7, there are 19) S_BUT run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 19 individuals (38 observations).

# 2*number of individuals in data set
N=38

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind38_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind38")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "S_BUT"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(!seg_sites == 0) %>%
  filter(pop == grp) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.15249, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.17 S_FRH

The maximum number of segregating sites per locus in the data set is 7, there are 7) S_FRH run individuals in the data set.

⊕Simulate 1000 independent replicates (loci) for 1 - 7 segregating sites under neutral model for 7 individuals (14 observations).

# 2*number of individuals in data set
N=14

# number of loci to simulate
REP=1000

# simulate neutral loci
for i in {1..7}
do

  scr/MS/ms $N $REP -s $i | scr/MS/sample_stats > results/Ind14_Sites${i}_neut_Rep1000.stats

done

The sample_stat script distributed along side MS calculates summary statistics including Tajima’s \(D\) for simulated set of loci.

⊕Generate 100 simulated data sets based on empirical data set, by determining the number of loci with 1 … x segregating sites and randomly drawing that number of loci from simulated data set (1,000 simulated loci per bin; draw with replacement) to generate null distribution consisting of the same number of loci with the same distribution of segregating sites per locus as the empirical data set.

# create a vector of stats files based on simulations
filenames <- list.files(path = "results/", pattern = "Ind14")

# create an empty list that will serve as a container to receive the incoming files
sim <-list()

# create a loop to read in your data
for (i in 1:length(filenames)){

  sim[[i]] <- read_delim(file.path("results", filenames[i]), delim = "\t",
                         col_names = c("t1", "pi",
                                       "t2", "seg_sites",
                                       "t3", "tajima_d",
                                       "t4", "thetaH",
                                       "t5", "H")) %>%
    select(pi, seg_sites, tajima_d, thetaH, H)

}

# add names
names(sim) <- filenames

# combine into dataframe
sim <- ldply(sim, data.frame) %>%
  rename(simulation = `.id`) %>%
  separate(simulation, into = c("temp", "t1"), sep = -6) %>%
  separate(temp, into = c("N_IND", "t2", "scenario", "N_REPS"), sep = "_") %>%
  mutate(N_REPS = as.numeric(str_extract(N_REPS, "[[:digit:]]+")),
         N_IND = as.numeric(str_extract(N_IND, "[[:digit:]]+"))) %>%
  select(-t1, -t2)

reps <- 100

# group to be compared
grp <- "S_FRH"

# get counts for number of segregating sites
emp <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites) %>%
  filter(pop == grp) %>%
  filter(!seg_sites == 0) %>%
  count(seg_sites)

# create null distributions
null_distributions <- list()

for (s in emp$seg_sites) {

      # determine number of loci with s segregating  sites
      t <- emp %>%
        filter(seg_sites == s)

      n_loci <- t$n

     # sample n_loci simulated loci (with replacement)
     null_distributions[[s]] <- sim %>%
        filter(seg_sites == s) %>%
        sample_n(size = n_loci, replace = TRUE)

}

null_distributions <- ldply(null_distributions, data.frame) %>%
    mutate(pop = grp)

taj <- read_delim("results/trib_run.hapdiv", delim = "\t") %>%
  select(locus, pop, seg_sites, tajima_d) %>%
  filter(!seg_sites == 0,
         pop == grp) %>%
  bind_rows(null_distributions) %>%
  mutate(scenario = ifelse(is.na(scenario), "empirical", scenario)) %>%
  select(-N_IND, -N_REPS, -pi, -thetaH, -H)

taj_all[[grp]] <- taj

Use K-S test statistic to identify simulated null distribution most similar to empirical data set (values closer to 0 indicate higher probability to have been draw from the same distribution).

Cumulative distribution curves of empirical data set (black) to simulated data set (grey). The black dashed lines indicates expected value for Tajima’s D under equilibrium assuming a beta distribution (the median should intersect at 0).

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.30666, p-value < 0.00000000000000022
## alternative hypothesis: two-sided

3.8.2.18 Comparison across all trib/run populations

Let’s compare the null and empirical distributions across all data sets.

Comparison of neutral and empirical distributions of Tajima’s D for each run in each tributary.

Table S10: Mean, median, 25% and 75% quantiles Tajima’s D forempirical and simulated data sets.

Difference in mean Tajima’s D between empirical and simulateddata sets.

Difference in mean Tajima’s D between empirical and simulateddata sets.

Typically the mean is lower than the median but we are more interested in genome-wide patterns so the median is more informative.

The median lower in the empirical data sets compared to the simulated data set under mutation drift equilibrium in the following populations.